home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / src / balance.cc < prev    next >
C/C++ Source or Header  |  1996-11-03  |  6KB  |  287 lines
  1. /*
  2.  
  3. Copyright (C) 1996 John W. Eaton
  4.  
  5. This file is part of Octave.
  6.  
  7. Octave is free software; you can redistribute it and/or modify it
  8. under the terms of the GNU General Public License as published by the
  9. Free Software Foundation; either version 2, or (at your option) any
  10. later version.
  11.  
  12. Octave is distributed in the hope that it will be useful, but WITHOUT
  13. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  14. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  15. for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with Octave; see the file COPYING.  If not, write to the Free
  19. Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  20.  
  21. */
  22.  
  23. // Written by A. S. Hodel <scotte@eng.auburn.edu>
  24.  
  25. #ifdef HAVE_CONFIG_H
  26. #include <config.h>
  27. #endif
  28.  
  29. #include <string>
  30.  
  31. #include "CmplxAEPBAL.h"
  32. #include "CmplxAEPBAL.h"
  33. #include "dbleAEPBAL.h"
  34. #include "dbleAEPBAL.h"
  35. #include "dbleGEPBAL.h"
  36.  
  37. #include "defun-dld.h"
  38. #include "error.h"
  39. #include "gripes.h"
  40. #include "help.h"
  41. #include "oct-obj.h"
  42. #include "utils.h"
  43.  
  44. DEFUN_DLD (balance, args, nargout,
  45.   "AA = balance (A [, OPT]) or [[DD,] AA] =  balance (A [, OPT])\n\
  46. \n\
  47. generalized eigenvalue problem:\n\
  48. \n\
  49.   [cc, dd, aa, bb] = balance (a, b [, opt])\n\
  50. \n\
  51. where OPT is an optional single character argument as follows: \n\
  52. \n\
  53.   N: no balancing; arguments copied, transformation(s) set to identity\n\
  54.   P: permute argument(s) to isolate eigenvalues where possible\n\
  55.   S: scale to improve accuracy of computed eigenvalues\n\
  56.   B: (default) permute and scale, in that order.  Rows/columns\n\
  57.      of a (and b) that are isolated by permutation are not scaled\n\
  58. \n\
  59. [DD, AA] = balance (A, OPT) returns aa = dd*a*dd,\n\
  60. \n\
  61. [CC, DD, AA, BB] = balance (A, B, OPT) returns AA (BB) = CC*A*DD (CC*B*DD)")
  62. {
  63.   octave_value_list retval;
  64.  
  65.   int nargin = args.length ();
  66.  
  67.   if (nargin < 1 || nargin > 3 || nargout < 0 || nargout > 4)
  68.     {
  69.       print_usage ("balance");
  70.       return retval;
  71.     }
  72.  
  73.   string bal_job;
  74.   int my_nargin;        // # args w/o optional string arg
  75.  
  76.   // Determine if balancing option is listed.  Set my_nargin to the
  77.   // number of matrix inputs.
  78.  
  79.   if (args(nargin-1).is_string ())
  80.     {
  81.       bal_job = args(nargin-1).string_value ();
  82.       my_nargin = nargin-1;
  83.     }
  84.   else
  85.     {
  86.       bal_job = "B";
  87.       my_nargin = nargin;
  88.     }
  89.  
  90.   octave_value arg_a = args(0);
  91.  
  92.   int a_nr = arg_a.rows ();
  93.   int a_nc = arg_a.columns ();
  94.  
  95.   // Check argument 1 dimensions.
  96.  
  97.   int arg_is_empty = empty_arg ("balance", a_nr, a_nc);
  98.  
  99.   if (arg_is_empty < 0)
  100.     return retval;
  101.   if (arg_is_empty > 0)
  102.     return octave_value_list (2, Matrix ());
  103.  
  104.   if (a_nr != a_nc)
  105.     {
  106.       gripe_square_matrix_required ("balance");
  107.       return retval;
  108.     }
  109.  
  110.   // Extract argument 1 parameter for both AEP and GEP.
  111.  
  112.   Matrix aa;
  113.   ComplexMatrix caa;
  114.   if (arg_a.is_complex_type ())
  115.     caa = arg_a.complex_matrix_value ();
  116.   else
  117.     aa = arg_a.matrix_value ();
  118.  
  119.   if (error_state)
  120.     return retval;
  121.  
  122.   // Treat AEP/GEP cases.
  123.  
  124.   switch (my_nargin)
  125.     {
  126.     case 1:
  127.       
  128.       // Algebraic eigenvalue problem.
  129.  
  130.       if (arg_a.is_complex_type ())
  131.     {
  132.       ComplexAEPBALANCE result (caa, bal_job);
  133.  
  134.       if (nargout == 0 || nargout == 1)
  135.         retval(0) = result.balanced_matrix ();
  136.       else
  137.         {
  138.           retval(1) = result.balanced_matrix ();
  139.           retval(0) = result.balancing_matrix ();
  140.         }
  141.     }
  142.       else
  143.     {
  144.       AEPBALANCE result (aa, bal_job);
  145.  
  146.       if (nargout == 0 || nargout == 1)
  147.         retval(0) = result.balanced_matrix ();
  148.       else
  149.         {
  150.           retval(1) = result.balanced_matrix ();
  151.           retval(0) = result.balancing_matrix ();
  152.         }
  153.     }
  154.       break;
  155.  
  156.     case 2:
  157.       {
  158.     // Generalized eigenvalue problem.
  159.  
  160.     // 1st we have to check argument 2 dimensions and type...
  161.  
  162.     octave_value arg_b = args(1);
  163.  
  164.     int b_nr = arg_b.rows ();
  165.     int b_nc = arg_b.columns ();
  166.       
  167.     // Check argument 2 dimensions -- must match arg 1.
  168.  
  169.     if (b_nr != b_nc || b_nr != a_nr)
  170.       {
  171.         gripe_nonconformant ();
  172.         return retval;
  173.       }
  174.       
  175.     // Now, extract the second matrix...
  176.     // Extract argument 1 parameter for both AEP and GEP.
  177.  
  178.     Matrix bb;
  179.     ComplexMatrix cbb;
  180.     if (arg_b.is_complex_type ())
  181.       cbb = arg_b.complex_matrix_value ();
  182.     else
  183.       bb = arg_b.matrix_value ();
  184.  
  185.     if (error_state)
  186.       return retval;
  187.  
  188.     // Both matrices loaded, now let's check what kind of arithmetic:
  189.  
  190.     if (arg_a.is_complex_type () || arg_b.is_complex_type ())
  191.       {
  192.         if (arg_a.is_real_type ())
  193.           caa = aa;
  194.  
  195.         if (arg_b.is_real_type ())
  196.           cbb = bb;
  197.  
  198.         // Compute magnitudes of elements for balancing purposes.
  199.         // Surely there's a function I can call someplace!
  200.  
  201.         for (int i = 0; i < a_nr; i++)
  202.           for (int j = 0; j < a_nc; j++)
  203.         {
  204.           aa (i, j) = abs (caa (i, j));
  205.           bb (i, j) = abs (cbb (i, j));
  206.         }
  207.       }
  208.  
  209.     GEPBALANCE result (aa, bb, bal_job);
  210.  
  211.     if (arg_a.is_complex_type () || arg_b.is_complex_type ())
  212.       {
  213.         caa = result.left_balancing_matrix () * caa
  214.           * result.right_balancing_matrix ();
  215.  
  216.         cbb = result.left_balancing_matrix () * cbb
  217.           * result.right_balancing_matrix ();
  218.  
  219.         switch (nargout)
  220.           {
  221.           case 0:
  222.           case 1:
  223.         warning ("balance: should use two output arguments");
  224.         retval(0) = caa;
  225.         break;
  226.  
  227.           case 2:
  228.         retval(1) = cbb;
  229.         retval(0) = caa;
  230.         break;
  231.  
  232.           case 4:
  233.         retval(3) = cbb;
  234.         retval(2) = caa;
  235.         retval(1) = result.right_balancing_matrix ();
  236.         retval(0) = result.left_balancing_matrix ();
  237.         break;
  238.  
  239.           default:
  240.         error ("balance: invalid number of output arguments");
  241.         break;
  242.           }
  243.       }
  244.     else
  245.       {
  246.         switch (nargout)
  247.           {
  248.           case 0:
  249.           case 1:
  250.         warning ("balance: should use two output arguments");
  251.         retval(0) = result.balanced_a_matrix ();
  252.         break;
  253.  
  254.           case 2:
  255.         retval(1) = result.balanced_b_matrix ();
  256.         retval(0) = result.balanced_a_matrix ();
  257.         break;
  258.  
  259.           case 4:
  260.         retval(3) = result.balanced_b_matrix ();
  261.         retval(2) = result.balanced_a_matrix ();
  262.         retval(1) = result.right_balancing_matrix ();
  263.         retval(0) = result.left_balancing_matrix ();
  264.         break;
  265.  
  266.           default:
  267.         error ("balance: invalid number of output arguments");
  268.         break;
  269.           }
  270.       }
  271.       }
  272.       break;
  273.  
  274.     default:
  275.       error ("balance requires one (AEP) or two (GEP) numeric arguments");
  276.       break;
  277.     }
  278.  
  279.   return retval;
  280. }
  281.  
  282. /*
  283. ;;; Local Variables: ***
  284. ;;; mode: C++ ***
  285. ;;; End: ***
  286. */
  287.